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1.  Introduction 

Intensive  efforts  are  underway  to  develop  small  PET  systems  dedicated  to  functional 
imaging  of  breast  cancer.  These  system  developments  will  still  require  sophisticated 
analyses  of  dynamic  PET  images  to  detect  small  lesions.  We  have  been  seeking  a 
clinically  practical  way  to  assist  conventional  visual  inspection  of  static  PET  images  with 
temporal  information  derived  from  a  dynamic  PET-FDG  imaging  sequence.  Modeling, 
extracting  and  exploiting  physiological  features  that  can  distinguish  normal  breast  tissues 
from  axillary  malignances  in  dynamic  PET  has  been  pursued  since  the  last  reporting 
period.  The  study  conducted  in  the  current  year  was  built  on  the  previous  year's  research, 
which  concentrated  on  the  development  and  testing  of  optimal  computer-aided  detection 
criteria  for  kinetic  feature  identification  from  equivocal  axillary  metastases  in  noisy  images. 
Rather  than  replace  conventional  visual  image  inspection,  our  goal  is  to  design  an 
intelligent  system  that  will  supplement  it. 

The  whole  proposed  study  consists  of  four  tasks:  Task  1:  Developing  the  mathematical 
formula  to  linearly  map  and  identify  the  physiological  features  contained  in  PET  dynamic 
sinogram  sequence  {Month  1-8),  Task  2:  Developing  the  schemes  for  objective  reduction 
of  dynamic  sinogram  data  guided  by  the  identified  TAG  subspaces  of  the  desired  signal 
(tumor)  and  the  interference  (normal  tissue  background  plus  noise)  {Month  4  -  12),  Task 
3:  Deriving  and  analyzing  statistical  hypothesis  test  criteria  to  test  the  presence  of  an 
axillary  metastasis  in  the  dynamic  images  reconstructed  from  the  compressed  sinogram 
data  {Month  13  -  24),  and  Task  4;  Clinical  Evaluation  {Month  13  -  36).  For  each  task, 
several  subtasks  were  defined  (see  the  SOW  in  the  grant  application  for  details). 

2.  Body 

During  this  annual  reporting  period,  our  efforts  were  mainly  focused  on  Tasks  3  and  4. 
These  two  tasks  were  those  originally  proposed  for  the  second  year  of  the  study.  After  last 
year’s  annual  report,  it  was  brought  to  our  attention  that  we  needed  to  obtain  approval 
from  the  U.S.  Army  Medical  Research  and  Material  Command  Institutional  Review  Board 
(US  Army  IRB).  Although  patients  whose  data  are  used  in  this  study  will  have  a  PET  scan 
regardless  of  their  participation  in  this  study,  they  still  must  be  prospectively  recruited  for 
“additional”  (temporal-based)  images.  Efforts  have  been  underway  for  the  last  6  to  8 
months  to  complete  the  respective  IRB  approvals  with  protocol/language  acceptable  to 
use’s  Institutional  Review  Board  and  the  US  Army  IRB.  Because  IRB  approval  is 
pending,  no  dynamic  imaging  data  have  been  acquired  nor  any  subjects  recruited  for  the 
extra  pictures  since  that  time.  Efforts  in  support  of  Task  4  have  been  limited  to  analyses  of 
dynamic  imaging  data  that  already  existed  in  the  patient  database  of  USC  PET  center. 
Further  work  on  Tasks  1  and  2  has  also  been  continued  through  this  year.  The  major 
accomplishments  of  this  activity  are  presented  as  follows.  The  technical  details  and 
related  results  can  be  found  in  Appendices. 

2.1  Evaluation  and  Validation  of  PET-FDG  tracer  kinetic  models  and  subspace 
representations  {Task  1  continuation). 
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Last  year  utilizing  a  literature  search,  we  reviewed  various  types  of  mathematical  models 
of  PET-FDG  kinetics.  We  selected  the  three-compartmental  (or  4-K  compartmental  model) 
to  use  in  these  analyses  [1]. 

Evaluation  and  validation  of  the  standard  PET-FDG  3-compartment  metabolic  models  in 
homogeneous  and  heterogeneous  tissues  was  performed  on  dynamic  PET  data  from  the 
image  database  at  the  USC  PET  center.  Among  them,  7  breast  cancer  and  5  lung  cancer 
cases  were  used  for  the  kinetic  analysis.  In  these  images  there  were  7  primary  breast 
tumors,  2  axillary  metastases,  5  primary  lung  tumors  as  well  as  2  lung  metastases.  All 
were  proven  by  pathology. 

Using  the  clinical  data,  the  hypothesis  that  the  kinetic  features  of  primary  tumors  and 
metastases  are  similar  but  they  are  different  in  malignant  and  normal  tissues  was  tested. 
For  each  patient  data,  ROIs  were  selected  in  well-defined  primary  tumor,  metastases  and 
normal  tissues,  respectively.  In  the  ROIs,  we  observed  that  the  radiotracer  uptake  within 
the  primary  tumors  or  metastases  increases  with  time,  while  the  TAG  decreases  for 
normal  tissues.  To  further  reveal  the  similarities  and  differences  in  kinetic  parameters,  two 
methods  were  employed  to  fit  the  time  activity  curves  measured  in  different  tissue  types 
with  the  three  compartmental  kinetic  models.  The  first  assumed  homogeneity  within  the 
ROIs  and  the  other  method  accounted  for  potential  heterogeneity.  In  the  homogeneous 
PET-FDG  modeling,  the  nonlinear  parameter  estimation  with  the  Newton-Raphson 
algorithm  was  used  to  estimate  the  rates  of  transport  of  FDG  in  plasma  to  tissue,  Ki,  k2, 
ks,  and  k4,  while  fitting  the  time  activity  curves.  The  corresponding  macro-parameters, 
which  define  the  physiological  factors  uniquely  associated  to  tissue  types,  were  computed 
from  the  estimations  of  k/s.  In  heterogeneous  tissue,  the  time  activity  curve  generated  by 
the  ROI  represents  the  mass-weighted  average  concentration  of  radioactivity  in  several 
different  "homogeneous"  tissues.  Thus,  the  physiological  factors  demonstrated  in 
heterogeneous  tissues  really  represent  those  from  several  homogeneous  tissues  [1,2]. 
The  multidimensional  maximum  likelihood  estimation  (MLE)  method  [7,8]  was  utilized  to 
estimate  the  macro-parameters  of  physiological  factors  in  heterogeneous  tissues.  We 
found  that  no  significant  "distance"  between  macro-parameter  clusters  of  primary  and 
metastases.  But  the  "distance"  of  macro-parameter  clusters  between  malignancy  and 
normal  tissue  was  significantly  different  (  >  0.006  and  0.001)  for  breast  and  lung  cancers, 
respectively.  We  also  found  that  the  use  of  heterogeneous  model  and  the  MLE  results  in 
improved  lesion  detectability,  specifically  with  respect  to  false  positive  and  false  negative 
findings,  although  the  mean  square  errors  generated  by  fitting  the  time  activity  curve  were 
about  the  same  for  the  two  models.  Since  the  initial  dataset  studied  is  too  small  to  make  a 
conclusive  decision,  the  assessment  of  the  two  models  and  their  associated  kinetic 
parameter  estimations  will  be  continued  in  the  third  year  when  more  clinical  data  will  be 
available. 

2.2  Identification  of  Physiological  subspace  directly  from  PET  projection  data 

(Continuation  for  Task  2) 

A  data  processing  scheme  was  developed  to  directly  identify  the  kinetic  parameters  and 
the  corresponding  subspace  features  of  the  lesions.  Instead  of  reconstructing  all  frames  of 
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dynamic  images  for  kinetic  analysis,  our  method  selectively  reconstructs  only  one  frame. 
This  frame  is  most  often  the  one  with  the  longest  acquisition  time  period  or  highest  signal- 
to-noise  ratio.  From  the  selected  image,  the  ROIs  can  be  defined,  which  provides  in  turn 
the  ROIs  for  the  sinogram  domain.  Different  than  the  image  data,  each  observation  in 
sinogram  contains  contributions  from  all  image  pixels  on  a  ray  line.  Thus,  filtering  has  to 
be  performed  in  order  to  remove  the  superimposed  contributions  from  other  pixels  before 
proceeding  to  estimation  of  subspace  parameters  in  the  sinogram  ROIs.  The  filtered  data 
provide  a  more  accurate  parameter  estimation.  A  filter  was  designed  which  turned  out  to 
be  a  local  filtered  backprojection  operator.  The  implementation  of  the  filter  in  C++  was 
completed. 

2.3  Performance  test  of  data  compression  algorithms  for  fast  dynamic  PET  image 
reconstruction  (Task  2  continuation) 

Three  rank-reduction  criteria  that  lower  the  dimensionality  of  dynamic  sinogram  in  time 
domain  and  speed  up  dynamic  image  reconstruction,  developed  in  the  first  year,  were 
further  evaluated  with  computer  generated  phantom  dynamic  data.  The  time  activity 
curves  of  the  primary,  metastatic  lesions  and  normal  tissues  in  the  phantom  data  were 
true  clinical  observations.  A  ROC  study  on  the  phantom  data  is  ongoing  to  assess  these 
three  algorithms  in  terms  of  lesion  detectability  (false  positive  and  false  negative  rates).  To 
date,  our  results  show  the  proposed  maximum  SNR  criterion  outperforms  the  conventional 
SVD  and  maximum  signal  energy  methods  in  terms  of  enhancing  lesion  to  background 
contrast  in  reconstructed  images.  This  is  because  that  the  maximum  SNR  method  best 
maintains  the  characteristics  of  the  lesion  time  activity  curve  after  compression,  while  both 
the  SVD  and  the  maximum  signal  energy  methods  [3,  9]  change  the  characteristics  of  time 
activity  curve  in  the  lesion  substantially. 

The  performance  assessment  of  the  compression  algorithm  with  clinical  data  will  continue 
to  be  pursued  in  the  third  year.  The  data  processing  scheme  and  the  related  software, 
described  in  section  2.2,  is  to  be  applied  to  estimate  kinetic  features  directly  from  PET 
sinogram.  The  estimated  features  will  be  used  in  the  developed  maximum  SNR  criterion  of 
data  compression. 

2.4  Development  of  statistical  hypothesis  test  criteria  for  the  computer-aided 
metastasis  detection  (Task  3) 

Using  the  subspace  kinetic  features  identified  from  well-defined  malignances  for 
computer-aided  detection  of  non-palpable  metastases  was  accomplished  by  the 
development  of  statistical  hypothesis  test  criteria  and  decision  rules.  These  criteria  are 
aimed  at  distinguishing  the  subspace  features  intrinsic  to  both  malignancy  and  normal 
tissues. 

The  generalized  likelihood  ratio  test  (GLRT)  is  a  standard  procedure  for  solving  statistical 
detection  problems  [5].  During  this  year,  three  types  of  GLRTs  were  derived  for  the 
different  assumptions  on  pixel  spatial  correlation.  The  first  GLRT,  that  is  an  extension  of 
the  algorithm  first  developed  in  [5],  assumes  that  pixels  in  an  arbitrarily  selected  ROI  are 
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spatially  uncorrelated.  Two  ditferent  hypothesized  data  models  were  used  to  characterize 
the  observations  of  large  lesions  and  small  lesions.  In  practice,  image  pixel  values  are  not 
statistically  uncorrelated  because  each  count  contributes  to  all  pixel  values  when  an 
image  is  reconstructed,  so  the  second  GLRT  adds  a  pre-whiting  procedure  to  the  first.  The 
spatial  covariance  matrix  required  for  the  pre-whitening  process  is  to  be  computed  using 
the  methods  given  in  [4,6].  Finally,  the  third  GLRT  extended  the  first  by  assuming  that  the 
spatial  inter-pixel  correlation  in  each  frame  has  the  same  structure  but  different  energy 
level. 

The  technical  details  in  sections  2.1  -  2.4  can  be  found  in  Dr.  C.C.  Huang’s  thesis,  which 
is  attached  to  this  report  as  Appendix  A. 

2.5  Comparison  and  evaluation  of  three  developed  hypothesis  test  criteria  (Tasks  3 
and  4) 

We  compared  three  developed  detection  criteria  using  three  different  methods  of  TAG 
evaluation  by  computer  generated  phantom  data.  The  first  method  studied  was 
conventional  averaging  over  an  ROI,  which  assumes  no  spatial  correlation  between 
pixels.  The  second  method  takes  the  spatial  inter-pixel  correlation  into  account.  The 
covariance  matrix  is  computed  frame  by  frame,  based  on  Huesman's  analytical  formulas 
and  Carson's  simplifications.  In  this  method,  the  TACs  for  individual  pixels  in  an  ROI  are 
first  spatially  decorrelated  by  the  covariance  matrix,  and  then  averaged  for  kinetic 
parameter  estimation.  The  third  method  assumes  that  the  spatial  inter-pixel  correlation  in 
the  /-th  frame  has  the  form  aiR,  i.e.,  each  frame  has  the  same  spatial  correlation  structure 
R,  but  a  different  energy  lever  cr,.  Thus,  the  covariance  matrix  computed  in  the  least  noisy 
frames  can  be  used  as  matrix  R  for  noisier  frames.  For  receiver-operator  curve  (ROC) 
study,  the  kinetic  features  estimated  from  these  three  methods  are  incorporated  into  the 
three  GLRTs  for  lesion  identification,  developed  from  the  generalized  likelihood  ratio 
principle  as  described  above.  The  rates  of  true  and  false  positives  are  counted  to  form 
ROC  curves.  We  used  50  sets  of  dynamic  phantom  simulations  with/without  artificially 
inserted  lesions.  These  lesions  were  assigned  clinically  measured  TACs  of  malignant  lung 
and  breast  tissues,  respectively.  The  results  show  that  Method  3  increased  specificity  and 
accuracy  by  20  and  11%,  respectively,  as  compared  to  the  other  two  methods.  No 
significant  difference  was  observed  between  Methods  1  and  2,  which  could  be  due  to 
inaccurate  covariance  matrix  estimation  introduced  by  the  simplification  made  in  the  early 
frames.  Further  study  of  this  will  continue  through  the  next  year. 

This  work  has  been  presented  in  SNM  annual  conference  2001,  Toronto,  Ca.,  June  2001, 
see  Appendix  D  for  details. 


3.  Key  Research  Accomplishments 

The  main  accomplishments  in  Year  2  are 

1.  Evaluation  of  PET-FDG  three  compartmental  models  in  homogeneous  and 
heterogeneous  tissues  using  clinical  breast  and  lung  cancer  patient  data,  validation 
of  the  subspace  representations  for  time  activity  curves  (TAC)  in  different  tissues. 
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such  as  primary  tumors,  metastases  and  normal  tissues,  etc.  and  quantification  of 
the  different  kinetics  between  normal  and  malignant  tissues  as  well  as  the  similar 
kinetics  between  primary  and  metastatic  lesion  kinetics. 

2.  Establishment  of  a  data  processing  scheme  to  estimate  kinetic  features  of  three 
compartmental  models  directly  from  PET  sinogram  and  implementation  of  the 
scheme  with  C++. 

3.  Design  of  a  filter  to  remove  the  contributions  of  unwanted  image  pixels  from 
sinogram  data  for  more  accurate  estimation  of  kinetic  features  from  PET  project 
data,  test  of  filter  performance  in  sinogram  data  and  extension  of  the  image-pixel- 
based  estimation  algorithms  to  a  scan-projection-based  approach. 

4.  Assessment  of  three  data  compression  algorithms,  developed  in  the  last  reporting 
period,  for  fast  dynamic  image  reconstruction  using  computer  generated  phantom 
data.  Comparison  of  the  receiver-operating  characteristic  (ROC)  curves  between 
the  original  and  compressed  FBP  dynamic  images. 

5.  Development  of  statistical  hypothesis  test  criteria  for  the  computer-aided 
metastasis  detection.  Derivation  of  three  statistical  hypothesis  test  criteria  by 
generalized  maximum  likelihood  ratio  principle  that  are  able  to  distinguish  the 
physiological  factors  (i.e.  subspace  features)  associated  to  malignances  and 
normal  tissues. 

6.  Assessment  of  the  three  different  hypothesis  test  criteria  by  performing  a  receiver 
operating  characteristic  (ROC)  study  on  computer  simulated  phantom  and  small 
size  of  clinical  data. 

7.  Accomplishment  of  the  respective  IRB  approvals  with  protocol/language  acceptable 
to  use’s  Institutional  Review  Board  and  the  US  Army  IRB. 

4.  List  of  Reportable  Outcomes: 

4.1  Publications 

1.  C.  C.  Huang,  "Computer-Aided  Lesion  Detection  in  Positron  Emission 
Tomography:  A  Signal  Subspace  Fitting  Approach",  Ph.  D.  Thesis,  submitted  to 
Electrical  Engineering  Department  of  USC,  May  2001. 

2.  X.  Yu,  C.  C.  Huang  and  P.  S.  Conti,  "Assessment  of  ROI-based  Time  Activity 
Analyses  in  Dynamic  PET  For  Oncology  "  SNM  Annual  Conference  2001,  Toronto, 
Ca.,  June  2001. 

3.  X.  Yu,  C.  C.  Huang  and  R.  Leahy,  "Applied  Subspace  Identification,  Filtering  and 
Identification  to  Lesion  Detection  in  Dynamic  PET  Images",  revising  for  a 
submission  to  IEEE  Trans,  on  Medical  Imaging,  December  2001. 
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4.2  Graduation 


Two  students  graduated  with  their  Ph.  D.  and  M.S.  degrees,  respectively.  Both  were 
partially  supported  by  this  award. 


5.  Conclusion 

In  Year  2,  study  of  physiological  feature  modeling,  tumor  identification  and  objective  rank- 
reduction  of  dynamic  PET  sinogram  for  fast  image  reconstruction  was  continued.  The 
mathematical  models  and  algorithms  developed  in  Year  1  were  validated  and  evaluated 
with  both  computer  generated  phantom  and  a  few  sets  of  clinical  data  selected  from  the 
patient  database  of  USC  PET  center.  We  found  that  there  were  no  significant  differences 
between  macro-parameters  of  primary  and  metastatic  tissue,  which  implies  there  are 
certain  kinetic  similarities  between  them.  In  contrast,  a  large  difference  was  noted 
between  malignant  and  normal  tissue  macro-parameters  for  both  breast  and  lung  cancers. 

The  major  focus  of  Year  2  was  to  develop  statistical  hypothesis  test  criteria  for  the 
computer-aided  detection  of  kinetic  features  in  metastases.  Three  types  of  generalized 
likelihood  ratio  tests  (GLRT)  were  derived  under  different  assumptions  on  pixel  spatial 
correlation.  We  compared  these  three  GLRTs  using  three  different  methods  of  TAG 
evaluation  and  computer  generated  phantom  data.  The  lesions  in  the  phantom  data  were 
equivocal  or  invisible  and  assigned  to  time  activity  curves  of  malignant  lung  and  breast 
tissues  clinically  measured  in  patient  data.  The  results  show  that  the  best  of  the  three 
increased  specificity  and  accuracy  by  20%  and  11%,  respectively,  compared  to  the  other 
two  methods. 

Results  of  Year  1  and  Year  2  study  indicate  that  it  is  potential  to  improve  lesion 
detectability  through  sophisticated  physiological  modeling  and  statistical  signal  processing 
techniques.  All  findings  in  theory  and  simulations  will  be  further  tested  and  evaluated  by  a 
larger  group  of  clinical  patient  data  with  proven  primary  breast  cancer  and  palpable 
axillary  metastases  in  Year  3. 
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Chapter  4 


Dynamic  Data  Formulation  and  Subspace 
Identification 

4.1  Overview 

In  computer-aided  lesion  detection  with  FDG-PET  dynamic  images,  the  TACs  de¬ 
fined  in  Chapter  3  are  useful  for  separating  lesions  from  normal  tissues  [29].  However, 
how  to  identify  the  exponential  subspace  functions  and  incorporate  the  identified 
subspaces  into  a  detector  remains  a  problem.  In  this  chapter,  based  on  the  dynamic 
temporal-spatial  data  structure,  signal-subspace-fitting  methods,  such  as  the  least 
squares  (LS)  method  and  the  singular  value  decomposition  (SVD)  method,  are  used 
to  estimate  parametric  (exponential  function)  and  nonparametric  (physiological  fac¬ 
tor)  subspaces,  respectively. 
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4.2  Dynamic  Data  Formulation 

One  of  the  advantages  in  the  FDG-PET  dynamic  study  is  that  the  dynamic  images 
can  be  formed  into  a  temporal-spatial  data  matrix  such  that  some  useful  signal  array 
processing  techniques  can  be  applied.  Let  y*  be  the  kinetic  data  measured  by  the 
PET  scanner  for  the  p-th  pixel  in  a  heterogeneous  ROI,  then 

Yk  =  Xfc  -f  rifc 

=  Eajt  nfc  (4-1) 


where  iip  denotes  the  noise  and  is  assumed  to  be  Gaussian.  Given  that  P  denotes 
the  number  of  pixels  in  the  ROI,  the  temporal-spatial  data  matrix  for  the  ROI  from 
a  PET  scanner  can  be  expressed  as 


Y  =  [xi,X2,---,xp] -1- [ni,n2,--*,np] 
=  EA-f-N 

where 


(4.2) 


A  =  [ai,a2,---,ap] 


N  =  [ni,n2,---,np]. 


(4.3) 
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In  array  signal  processing,  the  column  space  of  the  matrix  E  in  (4.2)  is  called 
the  signal  subspace  [69].  In  this  dissertation,  the  signal  subspaces  are  called  the 
lesion  and  normal  tissue  subspaces  which  are  denoted  as  (H)  and  (S),  respectively, 
where  (H)  is  the  subspace  spanned  by  the  columns  of  the  Nxp  matrix  H,  (S)  is 
the  subspace  spanned  by  the  columns  of  the  Nxt  matrix  S,  and  p  +  t  <  N.  Based 
on  the  compartmental  model  analysis,  the  subspaces  (H)  and  (S)  have  the  input 
plasma  subspace  in  common,  otherwise,  they  are  linearly  independent  in  which  case 
the  columns  of  the  concatenated  matrix  [HS]  span  the  (p  +  t)-dimensional  subspace 
(HS). 

4.3  Subspace  Identification  from  Know  Tissue 
Type 

Two  categories  of  subspace  identification  methods  are  used  in  this  dissertation, 
namely,  non-parametric  and  parametric  methods,  for  solving  the  signal-subspace- 
fitting  problem.  The  non-parametric  techniques  used  here  simply  decompose  the 
observed  data  to  extract  the  subspaces  using  a  singular  value  decomposition  method. 
The  parametric  method  is  used  to  estimate  the  parameters  of  interest  in  spectrum¬ 
like  functions,  e.g.,  least  squares  estimation  for  the  exponential  parameters  in  Eq. 
(3.7). 
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4.3.1  Non-Parametric:  Singulcir  Value  Decomposition 


The  non-parametric  methods  find  a  multidimensional  signal  subspace  by  a  simple 
eigen-decomposition  of  the  autocorrelation  matrix  or  a  singular  value  decomposition 
of  the  data  matrix.  The  SVD  method  is  used  in  this  dissertation  to  estimate  the 
physiological  factors  defined  in  Eq.  (3.17).  Given  that  the  coefficient  matrix  A 
in  Eq.  (4.2)  is  uncorrelated  with  the  noise  N  which  is  assumed  to  be  white,  the 
autocorrelation  matrix  of  the  data  Y  can  be  decomposed  as  follows: 

Ry  =  YY^ 

=  HAA^H  +  tr^I 

=  $.A,$r  +  (4.4) 

where  we  assume  that  AA^,  is  a  full  rank  matrix,  and  therefore  HAA^H^  can 
be  eigendecomposed  as  $.A$f  such  that  (H)=($.).  The  r  eigenvalues  of  the 
decomposition  combine  with  the  noise  covariance  to  form  the  r  x  r  diagonal  matrix 
As  =  A  +  cr^I,  with  the  eigenvalues  in  A,  arranged  in  decreasing  order.  The(Y- 
r)  X  (iV  -  r)  diagonal  matrix  An  contains  the  A  -  r  repeated  eigenvalues  al-  Thus 
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the  eigen-decomposition  of  the  data  autocorrelation  matrix  in  Eq.  (4.4)  results  in  a 
signal  subspace  and  a  noise-only  subspace 

4.3.2  Parametric:  Least  Squares  Estimation 

Finding  the  subspace  using  the  parametric  algorithm  is  equivalent  to  the  estimation 
of  the  exponential  basis  in  Eq.  (3.7),  which  is  similar  to  the  well-known  problem 
of  finding  the  direction  of  arrival  (DOA)  in  signal  subspace  processing  [69].  In  this 
dissertation,  the  least  squares  (LS)  method  is  applied  for  the  parametric  method 
and  the  best  fit  of  parameters  to  the  observed  data  are  chosen  to  constitute  the 
subspaces.  Subspace  fitting  using  the  LS  method  is  defined  as  [69]: 

H,  A  =  argmin||Y  —  HA]]!^  (4.5) 

where  ||  •  ||f  denotes  the  Frobenius  norm.  Because  the  subspace  fitting  problem 
is  separable  in  H  and  A  [69],  by  substituting  the  pseudoinverse  solution,  A  = 
(H^H)“^H^Y,  back  into  Eq.  (4.5),  the  following  equivalent  problem  is  derived: 

H  =  argmin||Y-H(H^H)-'H^Y|]^ 

=  argimn||P^Y||^ 

=  argmaxllPi/Yjl^ 

H. 
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(4.6) 


=  argmaxtr{PHYY^} 

H 

where  Ph  =  is  the  orthogonal  projection  matrix  that  projects  onto 

the  column  space  of  H  and  “tr”  means  the  trace  operator.  In  feature  identification 
using  parametric  algorithm  for  a  FDG-PET  dynamic  study,  the  only  unknowns  are 
the  exponential  parameters,  the  LS  method  becomes  to  search  the  parameters  over 
a  reasonable  range  for  the  tracer  kinetic  parameters. 

4.3.3  Subspace  Refining:  A  Subspace  Distance  Measure 

The  lesion  and  normal  tissue  subspaces  estimated  individually  by  the  parametric  or 
non-parametric  methods  from  the  known  types  of  tissue  data  generally  capture  the 
characteristics  of  lesion  and  normal  tissues,  respectively.  But  from  the  physiological 
compartmental  model  analysis,  the  resulting  two  subspaces,  containing  the  common 
input  plasma  subspace,  may  be  so  close  to  each  other  that  they  can  hardly  be 
separated.  Hence,  to  improve  the  lesion  detection  performance  (see  Chapter  5),  a 
subspace  correlation  method  is  used  to  refine  the  identified  subspaces. 

The  subspace  refining  procedure  selects  two  subsets  of  basis  vectors  (column 
vectors)  from  the  identified  subspaces  (H)  and  (S)  based  on  subspace  distance  max¬ 
imization  between  the  two  candidate  subsets,  subject  to  the  condition  that  the  LSE 
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are  less  than  given  values.  By  the  process  of  refining,  the  common  input  plasma 
function  will  be  excluded  from  both  lesion  and  normal  tissue  subspaces. 

By  a  subspace  distance  measure,  the  subspace  refining  can  be  described  as:  given 
the  identified  subspaces  (H)  and  (S), 


H',  S'  =  arg  max  distance  {H.S} 

H,S 


subject  to 


||Y-Piy/Yllp<ei 

||Y-Ps'YlU<e2 


(4.7) 


where  distance  {H,  S}  =  \Jl  -  rj  is  the  subspace  distance,  where  n  is  the  largest 
principal  correlation  coefficient  between  H  and  S  [21],  and  Ci  and  62  are  the  thresh¬ 
olds  set  for  achieving  the  fidelity  criterion. 


45 


Chapter  5 

Matched  Subspace  Detector 

5.1  Overview 

After  identifying  the  lesion  and  normal  tissue  subspaces  using  the  methods  described 
in  Chapter  4,  how.  to  incorporate  the  identified  subspaces  into  statistical  decision 
criteria  becomes  critical  and  will  be  addressed  in  this  Chapter.  The  generalized  like¬ 
lihood  ratio  test  (GLRT)  is  a  standard  procedure  for  solving  the  detection  problem 
[55,  75,  76,  77,  78].  The  GLRT  is  used  to  derive  the  matched  subspace  detector 
which  is  the  general  building  block  of  multi-rank  matched  filter  in  signal  processing. 
In  sonar  signal  processing,  the  matched  subspace  detector  is  also  called  a  matched 
field  detector  [54].  Generally,  the  matched  subspace  detectors  turn  out  to  be  a  ratio 
of  generalized  energy  detectors. 
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In  this  Chapter,  three  types  of  GLRTs  will  be  derived  under  different  assumption 
of  voxel  spatial  correlation.  The  first  GLRT  directly  applies  the  algorithm  devel¬ 
oped  in  [55]  by  assuming  that  pixels  are  spatially  uncorrelated.  Both  additive  and 
replacement  noise  models  are  used  to  characterize  the  observations  of  large  lesions 
and  small  lesions  in  this  approach.  Considering  that  image  pixel  values  are  not 
statistically  un correlated  because  each  count  contributes  to  all  pixel  values  [7],  the 
second  GLRT  adds  a  pre-whitening  procedure  to  the  first.  The  spatial  covariance 
matrix  required  for  the  pre-whitening  process  is  to  be  computed  using  the  methods 
described  in  Appendix  B  [7,  30,  31].  Finally,  the  third  detector  extends  the  GLRT 
criterion  to  a  multi-pixel  detection.  In  this  approach,  it  is  assumed  that  the  spa¬ 
tial  inter-pixel  correlation  in  each  frame  has  the  same  structure  but  different  energy 
level. 

5.2  Hypothesized  Data  Model 

In  this  dissertation,  depending  on  the  size  of  lesions,  two  data  models  are  applied 
for  the  detection  hypothesis,  namely,  the  replacement  model  and  the  superimposed 
model. 
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5.2.1  Replacement  Model 


Based  on  a  replacement  model,  which  is  used  to  model  a  large  lesion  in  FDG-PET 
images,  iwo  possible  hypotheses  regarding  how  the  PET  data  was  generated  are 
needed.  The  null  hypothesis  Hq  says  that  the  data  consists  of  a  sum  of  normal  tissue 
signal  xo  and  noise  no-  The  alternate  hypothesis  Hi  says  that  the  data  consist  of  a 
sum  of  lesion  signal  Xi  and  noise  ni.  That  is 

{  Ho  .  y  =  Xo  +  no  ^5 

\  If  1  :  y  =  Xi  +  ni 

where  the  signal  x,  is  assumed  to  obey  the  linear  subspace  model 

Xo  =  s<^,  (f>en\ 

Xi  =  H0,  H  €  deW’,  t<N-p.  (5.2) 

5.2.2  Superimposed  Model 

For  the  small  lesion  case,  partial  volume  effects  will  cause  the  lesion  signal  to  be 
mixed  with  the  background  normal  tissue  signal,  hence,  it  is  best  represented  by  a 
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superimposed  model.  To  model  the  small  lesion  case,  the  superimposed  model  is 
used  for  the  hypothesis  Hi  which  are  given  by 


Ho -.y  =  xoo  +  no 

Hi:y  =  Xi  +  xoi  +  ni 


(5.3) 


where  the  linear  model  for  the  Xi  (lesion  signal)  and  iij  are  the  same  as  those  defined 
in  the  replacement  model  and 


Xoo  =  S0O,  </>o  €  7^*,  under  hypothesis  Hq 

Xoi  =  under  hypothesis  Hi.  (5.4) 

5.3  GLRT:  Spatially  Uncorrelated  Assumption 

5.3.1  Replacement  Model 

Based  on  the  spatially  uncorrelated  assumption,  the  noise  n,  can  be  modeled  as 
normal  with  zero  mean  and  covariance  matrix  erf  I.  Then,  the  detection  problem  for 
the  replacement  model  becomes  a  test  of  the  distributions: 

i  Ho  :  y  =  Xo  +  no  :  A^[S^,  erfi] 

\  Hi  :  y^Ki+m:  Nine,  all]. 

The  likelihood  ratio  test  can  be  written  as 


/(y) 
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(5.6) 


Hence,  the  generalized  likelihood  ratio  test  (GLRT)  can  be  derived  by  substituting 


the  maximum  likelihood  estimate  (MLE)  of  n,-,  9  and  (f>: 


Ky)  = 


l{9,^l,y) 


llAill^  +  ^llnoll 


(5.7) 


The  MLE,  hi,  9,  and  can  be  written  as 


^  =  (S^S)-^S^y 

9  =  (H^H)-'H^y 

ho  =  y  -  =  y  -  S(S^S)”*S^y  =  P^y 

hi  =  y  -  =  y  -  H(H^H)-^H^y  =  P^y  (5-8) 


where  P^  is  the  orthogonal  projection  onto  the  complement  subspace  of  (H).  Be¬ 
cause  af  is  unknown  and  the  MLE  of  cr?  is  d?  =  it  is  more  convenient  to 

replace  the  GLRT  by  the  (iV/2)-root  GLRT.  Hence,  the  GLRT  for  the  replacement 
hypothesis  model  becomes 


7^r(y) 


=  [i{yr^  = 


Hfiolll 

Inilll 
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y^Psy 

y^P^y’ 


(5.9) 


♦ 


5.3.2  Superimposed  Model 

Under  the  spatially  uncorrelated  assumption,  the  detection  problem  for  superim¬ 


posed  model  is  a  test  of  the  distributions: 


Hq  :  y  —  Xqo  +  Ho  • 

Hi  :  y  =  Xi  -t-  Xoi  +  Hi  :  N[H.6  +  afl]. 


(5.10) 


Hence,  the  GLRT  can  be  expressed  as 


Ky)  = 


l{0Ji,ahy) 

K^Q,^l,y) 


|Ai1I^  +  2^I!AoI 


(5.11) 


The  MLE  ho  is  the  same  as  that  derived  in  the  replacement  model,  i.e.. 


ho  =  Psy- 


(5.12) 


For  the  hi,  because  the  MLE  of  6  and  (j)i  are  given  by 


=  [H  Sfy  =  ([H  S]^[H  S])"'  [H 


H^S  ■ 

-1 

S^H 

S^S 
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y 


(5.13) 


♦ 


(s’-pis)-'s^p^ 


and  from  Eq.  (A. 5)  in  Appendix  A,  then 


hi  =  y  —  —  S(^i 

=  y  -  H(H^P^H)-^H^Piy  -  S(S^P]^S)-^S^P^y 
=  y  —  Ej/sy  —  Esny 
=  PHsy 

where  Ehs  and  Ss//  are  called  the  oblique  projections  defined  m  Appendix  A.  There¬ 
fore,  the  ( A^/2)-root  GLRT  for  the  superimposed  model  becomes 

L.iy)  = 

_  (5.15) 

■  y^Pi^sy' 

Note  that  the  matrices  H  and  S  used  to  derived  the  above  GLRTs  for  both  re¬ 
placement  and  superimposed  models  are  assumed  independent  (after  subspace  angle 
refining).  In  Appendix  C,  we  showed  that  the  GLRTs  can  not  be  improved  by  fur¬ 
ther  projecting  the  lesion  subspace  (H)  onto  the  orthogonal  complement  subspace  of 
the  normal  tissue  subspace  (S).  The  GLRTs  for  the  superimposed  model  are  exactly 
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the  same  before  and  after  using  the  orthogonal  refined  subspaces.  For  the  replace¬ 
ment  model,  The  GLRT  is  degraded  after  using  the  refined  orthogonal  subspaces.  In 
summary,  the  GLRT  can  not  be  improved  if  two  orthogonal  subspaces  are  defined. 

5.4  Spatially  Correlated  Assumption 

5.4.1  Local  Statistics  Computation 

It  is  known  that  image  pixel  values  are  not  statistically  uncorrelated  because  each 
count  contributes  to  all  pixel  values  [7].  Hence,  a  decorrelation  for  the  pre-whitening 
process  is  important  before  applying  the  GLRT.  The  decorrelation  procedure  re¬ 
quires  the  computation  of  a  covariance  matrix  in  the  ROI.  The  covariance  matrix 
can  be  estimated  by  simply  averaging  the  samples  around  the  ROI.  However,  this 
sample-averaging  method  to  compute  the  local  statistics  can  be  poor  since  the  local 
statistics  assumption  required  may  be  invalid  [27,  79].  Calculation  accuracy  can 
be  guaranteed  by  using  the  formulas  developed  in  [7,  19,  30,  31]  which  exploit  the 
Poisson  data  model  for  a  PET  system.  The  linear  and  non-linear  covariance  com¬ 
putation  methods  for  the  FBP  and  MAP  reconstructed  images,  respectively,  are 
described  in  Appendix  B,  and  are  used  here  to  estimate  the  inter-pixel  covariance, 
frame  by  frame. 
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5.4.2  GLRT 


Because  of  the  spatially  correlated  assumption,  the  average  of  TACs  over  an  ROI 
may  enhance  rather  than  suppress  the  observation  noise.  Thus,  to  mitigate  this,  the 
TAC-averaging  is  preceded  by  a  spatial  decorrelation  operation  before  applying  the 
GLRT.  This  spatial  decorrelation  is  usually  called  a  pre-whitening  procedure. 

5.4.3  Multi-Pixel  GLRT 

In  this  subsection,  the  GLRT  criterion  is  extended  to  multi-pixel  GLRT  detection. 
In  this  approach,  it  is  assumed  that  the  spatial  inter-pixel  correlation  in  each  frame 
has  the  same  structure  but  different  energy  level.  Let  the  P-dimensional  vector  f,-  = 
[/,•(!),  /,(2),  •  •  • ,  /t( denote  the  spatial  pixel  data  for  the  i-th  frame  and  denote 
F  =  •  •  •  ,fiv].  We  assume  that  all  the  columns  of  the  matrix  F  are  independent. 

Hence  Y  =  [yi,y2>  •  ■  •  jyp]  =  derivation  of  the  multi-pixel  GLRT  for  the 

constrained  correlation  structure  is  based  on  a  replacement  model,  i.e., 

f  Ho  :  Y  =  Xo  -h  No  ^5 

\  Hi  :Y  =  Xi+Ni. 

The  GLRT  can  be  derived  as  (see  Appendix  B): 

_  l  +  (G-V^Rs)^P^(G-^/-Rs) 

-  1  +  (G-1/2Rs)TP^(G-1/2Rs) 
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where  G  =  RR^  —  (Rs)(Rs)^,  R  is  the  covariance  matrix,  and  s  =  (s^s)^/^s, 
s"^  =  The  matrix  M  is  the  correlation  structure  for  each  frame. 
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Chapter  6 

Experiment  and  Result 

In  this  Chapter,  experiments  and  results  for  dynamic  feature  identification  and  lesion 
detection  in  FDG-PET  dynamic  study  are  presented.  Both  clinical  and  phantom 
dynamic  data  will  be  used  in  this  dissertation.  For  the  clinical  dynamic  data,  two 
cases  will  be  demonstrated,  including  one  lung  cancer  dynamic  study  and  one  breast 
cancer  dynamic  study.  Each  study  contains  one  confirmed  primary  lesion  which 
will  be  used  for  lesion  subspace  estimation,  and  one  small  lesion  which  will  be  used 
for  the  subspace  fidelity  and  GLRT  demonstration.  For  the  breast  cancer  dynamic 
study,  there  is  one  “unknown”  small  lesion  which  can  not  be  confirmed  by  eye, 
but  is  successfully  enhanced  after  applying  the  GLRT.  For  the  phantom  data,  two 
phantoms,  one  with  five  artificial  lesions  inserted  and  the  other  with  only  pure  normal 
tissue  background  (no  lesion),  are  generated.  The  time  activity  curves  (TAG)  of  a 
lesion  and  normal  tissues  observed  from  the  clinical  lung  cancer  dynamic  data  are 
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used  to  simulate  the  dynamics  for  the  phantom.  To  calculate  the  local  statistics  for 
the  pre-whitening  purpose,  the  methods  [7,  19,  30,  31]  described  in  Appendix  B  will 
be  applied  for  covariance  computation.  A  receiver  operating  characteristic  (ROC) 
study  for  the  GLRT  performance  will  be  presented  based  on  the  dynamic  phantom 
study. 

6.1  Protocol  of  FDG-PET  Dynamic  Study 

The  clinical  FDG-PET  dynamic  data  was  acquired  with  a  Siemens/CTI  ECAT 
Model  953A  whole-body  PET  scanner.  This  device  provides  31  contiguous  transaxial 
image  planes  with  an  axial  field  of  view  of  10.8  cm.  The  nominal  intrinsic  resolution 
of  the  system  is  4  mm  in  all  3  dimensions.  Consecutive  detector  rings  are  separated 
by  tungsten  septa  in  order  to  reduce  scatter  noise.  Dynamic  data  was  acquired  from 
0  to  55  minute  post  injection.  The  final  10  minute  frame  can  be  used  in  place  of 
the  routine  static  study  used  in  purely  clinical  examinations.  The  dynamic  struc¬ 
ture  protocol  for  the  clinical  data  acquisition  is  listed  in  Table  6.1.  Based  on  this 
protocol,  an  FDG-PET  dynamic  study  contains  36  frames  for  a  fixed  plane  position. 
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Scan  Type 

Scan  Times 

Frame  Duration 

Dynamic  scan 

0-2  min 

15  sec/frame 

Dynamic  scan 

3-5  min 

30  sec/frame 

Dynamic  scan 

6-25  min 

1  min/frame 

Dynamic  scan 

26  -  45  min 

5  min/frame 

Static  scan 

45  -  55  min 

10  min 

Table  6.1:  Dynamic  data  acquisition  protocol 


6.2  Clinical  Study 


Two  clinical  FDG-PET  dynamic  studies,  one  lung  cancer  study  and  one  breast 
cancer  study,  will  be  used  in  this  dissertation.  Each  study  contains  one  confirmed 
primary  lesion  which  will  be  applied  for  lesion  subspace  estimation  and  one  small 
lesion  which  will  be  used  for  the  subspace  fidelity  and  GLRT  demonstration.  The 
last  frame  image  reconstructed  from  the  FBP  method  is  shown  in  Figure  6.1:  (a) 
and  (b)  for  the  lung  cancer  study  and  the  breast  cancer  image,  respectively. 

An  ROI,  indicated  by  a  rectangle,  can  be  drawn  over  a  segment  of  the  tissue 
on  the  dynamic  FDG-PET  images  and  then  duplicated  to  all  frames,  so  that  the 
activity  in  that  tissue  can  be  tracked  over  time.  The  advantage  of  the  ROI  method 
over  a  single  pixel  is  that  it  is  less  sensitive  to  noise  and  typically  results  in  a 
smaller  error  in  the  parameter  estimates.  Data  obtained  through  ROI  analysis  of 
dynamic  images  produces  a  tissue  time  activity  curve  (TAG).  This  curve  represents 
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the  counts/second/pixel  (or  count  s/second /ml,  if  calibration  data  is  available)  in  a 
given  region  as  a  function  of  time.  Figure  6.1:  (c)  and  (d)  shows  examples  of  TAG 
for  the  two  clinical  dynamic  studies. 


(c)  (d) 


Figure  6.1:  Clinical  dynamic  studies  used  in  this  dissertation:  (a)  the  lung  cancer 
image  (FBP),  (b)  the  breast  cancer  image  (FBP),  (c)  one  TAG  example  from  the 
lung  cancer  dynamic  study,  (d)  one  TAG  example  from  the  breast  cancer  dynamic 
study. 
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6.2.1  Lung  Cancer  FDG-PET  Dynamic  Study 

In  the  lung  cancer  FDG-PET  dynamic  study,  two  lesions,  one  primary  (large  and 
clearly  visualized)  and  the  other  a  metastasis  (smaller  and  barely  seen),  in  a  lung 
cancer  patient  with  multiple  metastases,  are  indicated  in  Figure  6.2  (a).  Three 
5x5  pixel  ROIs,  LI,  L2,  and  BG,  indicated  in  Figure  6.2  (a),  were  selected  from  the 
primary  lesion,  the  metastasis  lesion,  and  the  normal  lung  tissue,  respectively.  Their 
corresponding  TACs  are  plotted  in  Figure  6.2  (b)  where  the  TACs  from  the  LI  and 
L2  both  have  the  trend  to  increase  with  time  while  the  BG  TAG  does  not.  This 
demonstrates  the  property  of  FDG  kinetics  to  show  differential  uptake  in  lesions  and 

normal  tissue. 

6.2. 1.1  Subspace  Estimation 

The  TACs  observed  from  the  LI  and  BG  ROIs  in  Figure  6.2  (b)  were  used  to 
form  the  temporal-spatial  data  matrices  for  identifying  the  lesion  and  normal  tissue 
subspaces,  (H)  and  (S),  respectively. 

For  the  parametric  method,  the  LS  search  algorithm  was  used.  For  this  lung 
cancer  study,  the  TAG  subspaces  of  normal  lung  tissue  and  lung  cancer  estimated  by 
the  LS  method  were  found  to  be  spanned  by  the  exponential  basis  {e  } 

and  {1  — e"*®**,  respectively,  and  the  parameter  sets  were  {9i  =  0.00025, 
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TAGS  01  nOb  In  He  2Sn  tarns  o(  M  BP  Images 


Figure  6.2:  Test  image  reconstructed  from  FBP:  (a)  the  last  frame  of  the  dynamic 
FDG-PET  FBP  reconstructed  images,  where  three  ROIs  were  chosen:  LI  in  the  big 
lesion,  L2  in  the  small  lesion,  and  BG  in  the  background,  (b)  the  corresponding 
TACs  for  LI,  L2,  and  BG,  respectively. 

$2  =  0.005,  ^3  =  0.05}  and  {fSi  =  0.005,  ^2  =  0.0223,  /?3  =  0.00125}  for  BG  and  LI, 
respectively.  For  the  non-parametric  method,  we  applied  the  SVD  method,  where 
the  data  correlation  matrix  was  estimated  using  the  temporal-spatial  data  matrix 
formed  from  LI  and  BG  ROIs.  The  eigenvectors  corresponding  to  the  first  d\  and 
</2  significant  eigenvalues  of  the  two  estimated  correlation  matrices  were  chosen  to 
span  the  TAG  subspaces  of  the  lung  cancer  and  normal  lung  tissues,  respectively.  In 
our  study,  dl  =  d2  =  2. 

The  accuracy  (fidelity)  of  the  lesion  subspace  estimated  using  the  primary  lung 
tumor  was  examined  by  projecting  the  observed  LI  and  L2  TACs  onto  the  estimated 
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lesion  subspace.  Also,  the  accuracy  of  the  normal  tissue  subspace  estimated  using 
the  BG  ROI  was  examined  by  projecting  the  observed  BG  TAG  onto  the  estimated 
normal  tissue  subspace.  On  the  other  hand,  the  separability  of  the  subspace  was 
examined  by  projecting  the  Ll  and  L2  TACs  onto  the  normal  tissue  subspace,  and 
by  projecting  the  BG  TAG  onto  the  Ll  subspace.  The  results  from  the  LS  method 
and  the  SVD  method  are  shown  in  Figure  6.3.  It  was  demonstrated  that  the  TAG 
of  the  metastatic  lesion  can  be  represented  by  the  subspace  features  extracted  from 
the  big  lesion  with  a  high  accuracy  in  both  methods,  and  that  a  large  separability 
exists  between  the  lesion  subspace  and  the  normal  lung  tissue  subspace. 

6.2.1.2  GLRX 

A  test  region  which  includes  the  small  lesion  and  some  part  of  the  primary  lesion 
was  selected  as  indicated  by  a  rectangle  in  Figure  6.4  (a).  The  GLRT  (white  noise 
assumption)  with  either  a  superimposed  or  a  replacement  data  hypothesis  model 
in  Eqs.  (5.9)  and  (5.15)  was  used  to  perform  a  multiple  pixel  test  for  each  5x5 
pixel  ROI  contained  in  the  test  region.  The  corresponding  GLRT  result  is  shown  in 
Figure  6.4  (c)-(e)  which  indicates  that  the  matched  subspace  detection  incorporated 
with  the  lesion  and  normal  tissue  subspaces  extracted  by  the  SVD  method  and  the 
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large  TUMOR  LUNG  (BACKGROUND)  SMAa  TUMOR 


(d)  (e)  (f) 


Figure  6.3:  Comparison  of  accuracy  and  separability  of  extracted  subspaces  of  TACs 
in  lesions  and  normal  tissue  by  the  SVD  method:  (a)  TAG  in  LI  represented  by 
the  normal  tissue  subspace,  (b)  TAG  in  BG  represented  by  the  lesion  subspace, 
(c)  TAG  in  L2  represented  by  the  lesion  and  normal  tissue  subspaces,  respectively. 
Comparison  of  accuracy  and  separability  of  extracted  subspaces  of  TACs  in  lesions 
and  normal  tissue  by  the  LS  method:  (d)  TAG  in  LI  represented  by  the  normal 
tissue  subspace,  (e)  TAG  in  BG  represented  by  the  lesion  subspace,  (f)  TAG  in  L2 
represented  by  the  lesion  and  background  subspaces,  respectively. 


63 


LS  method  significantly  increases  the  lesion-to-normal  tissue  contrast  for  the  small 
lesion  relative  to  the  original  FBP  image  in  Figure  6.4  (b). 

Another  test  region  was  selected  to  include  the  heart  area  in  order  to  demonstrate 
suppression  of  activity  in  the  heart  by  the  GLRT.  Figure  6.5  (a)  shows  the  selected 
test  region  in  the  lung  cancer  dynamic  study,  where  the  primary  lesion,  the  small 
lesion,  and  the  heart  area  were  all  included.  Figure  6.5  (b)  shows  the  original  3-D 
mesh  of  the  test  region.  The  GLRT  of  the  test  region  is  shown  in  Figure  6.5  (c), 
where  the  heart  activity  is  suppressed  significantly  after  applying  the  GLRT. 

In  order  to  apply  the  whitening  process  for  the  spatially  uncorrelated  FBP  data, 
we  used  the  method  proposed  by  Carson  [7]  to  calculate  the  spatial  covariance  ma¬ 
trix.  Carson’s  approximation  formula  assumes  that  the  ROI  is  smooth  and  the  pixel 
variances  inside  the  ROI  are  similar.  Therefore,  the  covariance  matrix  is  just  a  func¬ 
tion  of  the  filter  and  the  location  of  the  pixels.  The  advantage  of  Carson’s  formula 
is  that  it  does  not  need  access  to  the  original  projection  data  and  is  convenient  for 
routine  use.  See  Appendix  B  for  more  details.  Figure  6.6  shows  the  results  of  the 
GLRT  for  the  FBP  image  after  applying  the  estimated  covariance  matrix. 
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3-0  MESH  d  TEST  REGION 


(c) 


(d) 


(e) 


Figure  6.4:  (a)  The  last  frame  of  the  dynamic  FDG-PET  FBP  reconstructed  images, 
where  a  rectangle  indicates  a  test  region  containing  a  part  of  the  big  lesion  and  the 
whole  small  lesion,  (b)  3-D  mesh  of  the  test  region  indicated  in  (a),  (c)  GLRT  with 
replacement  model  for  TAG  subspaces  obtained  by  the  SVD  method,  (d)  GLRT  with 
replacement  model  for  TAG  subspaces  obtained  by  the  LS  method,  (e)  GLRT  with 
superimposed  model  for  TAG  subspaces  obtained  by  the  LS  method. 
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Figure  6.5:  (a)  The  test  region  indicated  by  a  rectangle  in  the  lung  cancer  study, 
where  the  primary  lesion,  the  small  lesion,  and  the  heart  area  were  all  included,  (b) 
3-D  mesh  of  the  original  test  region,  (c)  The  result  of  the  GLRT. 


6.2. 1.3  Subspace  Refining  Method 

For  the  proposed  subspace  refining  algorithm  described  by  Eq.  (4.7),  we  estimated 
four  LS  TAG  bases  for  both  lesion  and  normal  tissue  subspaces,  and  used  the  sub¬ 
space  refining  algorithm  to  retain  only  three  bases  for  both  subspaces.  Figure  6.7 
(a)  and  (b)  shows  the  fidelity  between  the  estimated  TAG  by  the  identified  sub¬ 
spaces  and  the  observed  TAG  without  subspace  refining  (four  bases)  and  Figure  6.7 
(c)  shows  the  GLRT.  The  corresponding  results  after  subspace  refining  (three  bases 
retained)  are  shown  in  Figure  6.7  (d),  (e),  and  (f),  respectively,  where  not  only  the 
fidelity  was  preserved,  but  also  the  GLRT  performance  is  significantly  improved. 

For  the  orthogonal  projection,  we  already  showed,  see  Appendix  G,  that  the 
refined  subspaces  by  the  orthogonal  projection  can  not  improve  the  GLRT  with  the 
original  estimated  subspaces.  The  experimental  results  are  as  follows.  Figure  6.8 
(a)  and  (b)  show  the  test  region  chosen  and  its  corresponding  3-D  mesh.  Figure  6.8 
(c)  and  (d)  show  the  GLRTs  for  the  superimposed  model  before  and  after  using  the 
orthogonal  refined  subspaces,  where  the  two  GLRTs  were  exactly  the  same.  Then 
for  the  GLRT  of  the  replacement  model.  Figure  6.8  (f)  shows  that  the  GLRT  using 
the  refined  orthogonal  subspaces  had  poor  performance  compared  to  Figure  6.8  (e), 
where  the  original  subspaces  were  applied  for  the  GLRT. 
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Figure  6.7:  Subspace  refining  using  the  subspace  distance  measure.  Without  sub¬ 
space  refining:  (a)  lesion  TAC  and  (b)  normal  tissue  TAC  both  show  the  fidelity 
of  subspaces  estimated  from  observed  TAC,  (c)  GLRT.  With  subspace  refining:  (d) 
lesion  TAC  and  (e)  normal  tissue  TAC  shows  the  fidelity  of  subspace  refined.  It 
is  clear  that  the  refined  subspaces  maintain  the  fidelity  requirement,  (f)  GLRT 
performance  is  improved  by  using  the  refined  subspaces. 
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(e) 


(f) 


Figure  6.8:  Subspace  refining  using  the  orthogonal  projection,  (a)  test  region  inside 
the  lung  cancer  image,  (b)  original  3D  mesh.  GLRT  for  the  superimposed  model:  be¬ 
fore,  (c),  and  after,  (d),  the  orthogonal  subspace  refining.  GLRT  for  the  replacement 
model:  before,  (e),  and  after,  (f),  the  orthogonal  subspace  refining. 
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6,2.1,4  MAP  Reconstructed  Image 


We  also  tested  the  GLRT  for  the  lung  cancer  FGD-PET  dynamic  study  by  the  MAP 
reconstructed  algorithms  [42,  18,  43,  24].  For  the  non-white  noise  case,  the  fast 
covariance  calculation  method  proposed  by  Jinyi  and  Leahy  [31]  for  the  MAP  image 
was  used  here  for  the  decorrelation  purpose.  This  covariance  calculation  method 
was  summarized  in  Appendix  B.  We  used  the  same  subspaces  estimated  from  the 
FBP  image.  The  results  of  GLRT  are  shown  in  Figure  6.9. 

6.2.2  Breast  Cancer  FDG-PET  Dynamic  Study 

In  this  case,  we  examined  a  clinical,  36-frame,  dynamic  FDG-PET  study  of  a  patient 
with  a  primary  breast  cancer  and  a  known  axillary  metastasis.  The  primary  cancer 
was  clearly  visualized  in  the  left  breast,  while  the  smaller  axillary  metastasis  lymph 
node  was  marginally  visualized,  see  Figure  6.10  (a)-(b).  TACs  obtained  in  the  ROI 
placed  on  these  tumors  and  contralateral  normal  breast  tissue  are  shown  in  Figure 
6.10  (c).  The  TACs  showed  that  FDG  radiotracer  uptake  within  the  lesions  increases 
with  time,  while  the  TAG  decreases  for  normal  tissue. 

The  subspaces  for  the  breast  lesions  and  the  normal  tissues  were  identified  using 
the  SVD  method.  The  subspace  fidelity  and  separability  were  demonstrated  by  fit¬ 
ting  the  metastasis  lymph  node  ROI  onto  the  subspaces  estimated  from  the  primary 
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Figure  6.9:  The  MAP  reconstructed  lung  cancer  FDG-PET  image,  (a)  The  test 
region  indicated  by  a  rectangular  for  the  GLRT.  (b)  3-D  mesh  of  the  original  test 
region,  (c)  GLRT:  white  Gaussian  noise,  (d)  GLRT:  non-white  Gaussian  noise. 


cancer  ROI  and  the  normal  breast  tissue  ROI.  Figure  6.11  (d)  shows  that  the  sub¬ 
space  of  the  lesion  TACs  is  approximately  invariant  between  different  tumors  and 
the  subspaces  of  the  normal  tissue  and  the  lesion  TACs  were  different.  For  a  test 
region,  right  rectangle  in  Figure  6.11  (a),  containing  the  small  lesion,  applications  of 
the  GLRT  with  replacement  model  showed  that  the  lesion-to-normal  tissue  contrast 
was  enhanced  relative  to  the  original  FBP  image.  Figure  6.11  (f),  while  for  a  normal 
tissue  test  region,  left  rectangle  in  Figure  6.11  (a),  the  GLRT  shows  no  enhancement. 

An  interesting  GLRT  result  for  this  breast  cancer  dynamic  study  was  the  en¬ 
hancement  of  an  unknown  lesion  which  could  not  be  confirmed  by  eye.  This  un¬ 
known  lesion  was  roughly  located  (confirmed  by  a  later  scan)  inside  the  rectangular 
region  shown  in  Figure  6.12  (a)  and  (b),  but  was  undiscernible  to  human  eyes  even 
using  a  MAP  reconstruction  method.  This  patient  did  not  have  a  follow-up  scan  un¬ 
til  18  months  later.  From  the  newly  acquired  data,  the  “new”  lesion  was  confirmed. 
The  results  of  GLRT  (whit-noise  assumption)  based  on  the  earlier  FDG-PET  scan 
of  this  patient  showed  that  this  unknown  lesion  was  detected  successfully  in  both 
MAP  and  FBP  images,  see  Figure  6.12  (e)  and  (f),  when  compared  to  the  3-D  mesh 
of  the  test  region,  see  Figure  6.12  (c)  and  (d). 
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Figure  6.10:  Breast  cancer  study:  (a)  a  plane  containing  a  primary  cancer  indicated 
by  an  arrow,  (b)  a  plane  containing  a  metastatic  lymph  node  indicated  by  an  arrow, 
(c)  the  TACs  from  the  primary  cancer,  the  metastatic  lymph  node  and  normal  breast 
tissues  in  an  ROI  indicated  with  a  rectangle  in  (b). 

6.3  Phantom  Dynamic  Study 

6.3.1  Computer  Simulated  Phantom  Data 


Two  computer  simulated  phantoms,  one  with  five  artificial  lesions  at  the  known 
locations  and  the  other  with  no  lesions,  were  used  in  the  phantom  dynamic  study, 
see  Figure  6.13  (a)“(b).  The  TACs  of  the  normal  and  malignant  tissues  measured 
from  the  clinical  lung  cancer  FDG-PET  dynamic  study  were  applied  to  simulate  the 
lesion  and  normal  tissue  kinetics  behavior,  see  Figure  6.13  (c). 

During  the  forward  projection  of  the  phantom  dynamic  image,  the  total  count 
of  projection  data  in  each  frame  was  scaled  to  approach  the  total  count  of  the 
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Figure  6.11:  Breast  cancer  dynamic  study:  (a)  two  selected  test  regions:  test  region  1 
containing  the  normal  breast  tissue  and  test  region  2  containing  a  metastatic  lymph 
node  in  the  FBP  image  (same  as  (b)),  (b)-(c)  3-D  mesh  of  test  region  1  and  test 
region  2  indicated  in  (d).  (d)  TAG  of  the  metastatic  lymph  node  represented  by  the 
lesion  subspace  and  normal  breast  tissue  subspace  identified  by  the  SVD  method, 
respectively,  (e)-(f)  the  GLRT  results  on  test  region  1  and  2,  respectively. 
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(c)  (d) 


(e)  (f) 


Figure  6.12:  A  GLRT  result  of  enhancing  an  “unknown”  lesion  in  the  breast  cancer 
FDG-PET  dynamic  study.  The  test  region  containing  the  unknown  lesion  in  (a) 
FBP  and  (b)  MAP.  The  3-D  mesh  of  the  test  region:  (c)  FBP  and  (d)  MAP.  GLRT 
(white  noise):  (e)  FBP  and  (f)  MAP. 


(a)  (b)  (c) 


Figure  6.13:  Left:  lesion  phantom  image  with  5  lesions,  middle:  normal  tissue 
phantom  with  no  lesions,  right:  TACs  for  lesion  and  normal  tissues. 

corresponding  frame  in  a  clinical  dynamic  study.  This  will  be  called  the  low  count 

dynamic  study  in  this  research,  while  the  total  count  of  projection  data  for  a  high 

count  study  is  twice  of  the  total  count  in  a  clinical  dynamic  study.  During  the 

forward  process,  the  sinogram  data  were  corrupted  with  Poisson  noise  and  blurred 

kernels  to  simulate  the  noise  process  and  system  resolution.  FBP,  MAP,  and  OSEM 

reconstruction  methods  were  used  for  dynamic  phantom  image  reconstruction.  The 

images  (last  frame)  reconstructed  by  the  FBP,  MAP,  and  OSEM  methods  and  the 

corresponding  observed  TACs  of  normal  and  lesion  tissues  in  the  phantom  dynamic 

images  are  presented  in  Figure  6.14,  where  the  reconstructed  lesion  images  looked 

like  normal  tissue  images  and  the  artificial  lesions  were  hardly  confirmed. 
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The  GLRT  results  for  the  FBP  dynamic  images  are  shown  in  Figure  6.15  (high 
count  data)  and  6.16  (low  count  data),  and  the  GLRTs  for  the  MAP  low  count 
dynamic  phantom  images  were  shown  in  Figure  6.17.  The  GLRT  results  for  both 
high  count  data  and  low  count  data  improved  the  visualization  of  the  artificial  lesion 
locations. 

6.3.2  Receiver  Operating  Characteristic  (ROC)  Study 

An  ROC  study  was  applied  to  the  phantom  dynamic  study  for  comparing  the  perfor¬ 
mance  of  detectors.  Fifty  lesion  phantom  studies  (with  250  known  lesions)  and  fifty 
normal  tissue  phantom  dynamic  studies  were  conducted  in  a  Monte  Carlo  simula¬ 
tion.  The  rates  of  true  and  false  positive  decisions  made  from  the  three  GLRT  tests 
were  counted,  respectively.  The  ROC  curves  were  estimated  by  using  the  software 
package  ROCKIT  0.9B  developed  by  Professor  C.E.  Metz,  University  of  Chicago. 

The  resulting  ROC  curves  from  the  FBP  reconstructed  dynamic  phantom  images 
shows  that  the  GLRT  with  spatial  decorrelation  (with  area  under  curve  Az=0.92) 
are  close  to  that  for  the  GLRT  with  a  white-noise  (Az=0.93),  see  Figure  6.18  (a). 
The  ROC  curve  from  the  multi-pixel  GLRT  further  improved  the  Az  (=0.98). 
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(j)  (k)  (1) 

Figure  6.14:  Left  column:  lesion  phantom  image  with  5  lesions,  middle  column: 
normal  tissue  phantom  with  no  lesions,  right  column:  TACs  for  lesion  and  nor¬ 
mal  tissues.  From  top  row  to  bottom:  original  phantom,  FBP,  OSEM,  and  MAP 
reconstructed  images. 
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Figure  6.15:  The  GLRT  for 
white  noise  assumption,  (d] 


UMPpfwnkmM^)* 
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(b) 


(c) 


(d) 


Figure  6.17:  The  GLRT  for  an  MAP  phantom  image  (low  count  data):  (a)  a  test 
region  indicated  by  a  rectangular,  (b):  GLRT  (white  noise  assumption),  (c):  GLRT 
(non-white  noise  assumption),  (d):  multi-pixel  GLRT. 
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We  also  compared  the  GLRT  performance  for  different  reconstruction  algorithms. 
The  measured  Az  for  our  ROC  curves  can  be  viewed  as  a  computer  observers  for  dif¬ 
ferent  reconstruction  algorithms.  A  reconstruction  algorithm  can  be  claimed  “bet¬ 
ter”  if  it  achieves  a  higher  Az.  Figure  6.18  (b)  shows  that  the  Az  (MAP)  >  Az 
(OSEM)  >  Az  (FBP)  which  can  be  interpreted  as  the  iterative  reconstruction  meth¬ 
ods  (MAP  and  OSEM)  outperform  the  FBP  method  not  only  in  image  visualization, 
but  also  in  the  GLRT  detection  performance. 


(b) 


Figure  6.18:  ROC  curves  for  the  low  count  data  simulation:  (a)  FBP:  GLRT  (white 
noise,  non-white  noise)  vs  multi-pixel  GLRT,  where  Az  (multi-pixel  GLRT)  >  Az 
(GLRT  with  non-white  noise)  ~  Az  (GLRT  with  white  noise),  (b)  GLRT  with  white 
noise  assumption:  Az  (MAP)  >  Az  (OSEM)  >  Az  (FBP). 


(a) 
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Chapter  7 


Conclusion 

In  order  to  improve  early  detection  of  small  lesions  using  FDG-PET,  computer-aided 
lesion  detection  algorithms  have  been  developed  in  this  dissertation  to  assist  visual 
inspection  in  tumor  detection  using  dynamic  FDG-PET.  In  our  approach,  we  hy¬ 
pothesized  that  the  spatial  and  temporal  metabolic  features  available  from  dynamic 
FDG-PET  images  are  identifiable  in  a  clinically  practical  way  and  useful  for  lesion 
detection.  We  showed  that  malignancies  can  be  distinguished  from  normal  tissues 
on  the  basis  of  their  rates  of  FDG  uptake  in  terms  of  time  activity  curves,  given 
that  lesions  and  normal  tissues  often  differ  in  the  rate  of  radiotracer  accumulation 
or  disappearance.  Since  the  rates  of  FDG  uptake  in  tissues  are  often  characterized 
by  a  set  of  physiological  factors  and  the  time  activity  curves  can  be  represented  as 
a  linear  combination  of  the  physiological  factors,  we  used  the  subspace  spanned  by 
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the  physiological  factors  as  a  key  feature  to  distinguish  normal  tissues  from  malig¬ 
nancies.  By  forming  the  dynamic  FDG-PET  images  into  useful  spatial-temporal 
data  matrices,  we  were  able  to  estimate  the  tumor  and  normal  tissue  subspaces 
by  signal  subspace  fitting  methods  that  are  widely  used  in  array  signal  processing. 
Then,  applying  the  generalized  maximum  likelihood  ratio  principle,  we  adapted  and 
generalized  the  matched  subspace  detection  techniques  for  small  lesion  detection  in 
dynamic  FDG-PET.  We  demonstrated  that  with  the  matched  subspace  detectors, 
the  subspace  feature(s)  extracted  from  visible,  large  lesions  can  be  employed  to  con¬ 
firm  suspected,  but  not  unequivocally  identifiable,  small  lesions.  Results  from  both 
clinical  dynamic  FDG-PET  studies  of  patients  with  breast  or  lung  cancer  and  dy¬ 
namic  phantom  data  showed  that  the  physiological  subspaces  are  straightforward 
to  identify  and  distinct  for  known  tumors  and  normal  tissues,  and  are  valuable  for 
improving  lesion  diagnosis  when  combined  with  the  matched  subspace  detection  cri¬ 
teria.  Therefore,  the  detection  methods  proposed  in  this  dissertation  are  promising 
for  improving  early  detection  of  small  lesions.  Based  on  these  results,  it  is  our  be¬ 
lief  that  such  detection  methods  can  be  used  to  reduce  costly  multiple  diagnostic 
procedures  and  to  guide  surgical  intervention. 
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Appendix  A 

Oblique  Projection 


Let  (H)  denote  the  subspace  spanned  by  the  matrix  H.  If  the  subspaces  (H)  (p- 
dimension)  and  (S)  (t-dimension)  are  linearly  independent,  then  the  columns  of  the 
concatenated  matrix  [HS]  span  a  (p  +  t)-dimensional  subspace  (HS).  The  typical 
orthogonal  projection  of  a  TAG  y  (e  Tl^)  onto  (HS)  is  denoted  by  PhsY,  where 
Vhs  is  the  orthogonal  projection  in  the  subspace  (HS)  and  represented  as 

Phs  =  [HS]  ([HSj^iHS])"'  [HS]^.  (A.l) 


We  can  further  re-express  Phs  as 


Phs 


H  S 

■  H^H  H^S  ' 
S^H  S^S 

I _ 1 

H  0 

H^H  H^S 
S^H  S^S 

1 - 1 

S  H 

1 _ 1 

r  1 

H^H  H^S‘ 

-1 

H^ 

+ 

OS] 

S^H  S^S  _ 

.  s"’  . 
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(A.2) 


A 


Eh5  +  Esi/ 


where 


AT  T 

■  H^S  ' 

-1 

Ehs  =  [ho] 

S^H  S^S 

A  r 

H^H  H^S' 

-1 

Esh  =  [0  S] 

S^H  S^S 

The  two  projections,  E//5  and  Es^y,  are  called  the  oblique  projections  [6]  with  the 
property  that  Ehs  and  Es^  have  respective  range  spaces  (H)  and  (S)  and  respective 
null  spaces  (S)  and  (H),  i.e., 

EJ/5H  =  H,  EffsS  =  0, 

E5//S  =  S,  Esfl^H  =  0.  (A. 3) 


The  simplification  version  of  the  oblique  projections  E/rs  and  Esh  can  be  written 
as 


Esh  =  S(S^P^S)-^S"’P^.  (A.4) 
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Hence, 


Phs 


Ej/s  +  E5H 


H(H^PiH)-'H’'P^  +  S(S’'PiS)-'S’'P 


(A.5) 
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Appendix  B 

Local  Statistics  in  PET  Images 

B.l  Local  Statistics  Calculation  in  PET  FBP 
images 

The  analytic  computation  (AC)  for  the  local  statistics  of  a  PET  FBP  image  was 
originally  developed  by  Huesman  [17].  By  the  independence  of  observation  data  and 
the  linear  properties  of  FBP,  AC  estimates  the  (co)variance  of  PET  image  pixels 
from  the  projection  data  statistics  via  some  scaling  factors.  Huesman’s  algorithm 
calculates  local  statistics  directly  from  the  projection  data  without  image  recon¬ 
struction.  Recently,  Carson  et  al.  [7]  developed  an  approximation  formula  for  the 
variance  of  PET  ROI  values,  which  accounted  the  radioactivity  distribution,  atten 
nation,  random,  scatter,  deadtime,  detector  normalization,  scan  length,  decay,  and 
reconstruction  filter.  By  the  assumption  that  the  pixel  variances  of  any  two  pixels 
are  similar,  and  that  the  product  of  attenuation,  normalization,  wobbling  and  raw 
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data  is  relatively  uniform  for  those  projection  lines  that  substantially  contribute  to 
the  sum,  they  derived  a  simple  formula  for  the  correlation  coefficient  of  any  two 

pixels,  which  is  independent  of  the  raw  data. 

Let  the  raw  projection  data  in  PET  be  an  /-dimensional  column  vector  p  = 
[pf,  PnJ^5  where  I  is  the  product  of  the  total  number  of  angles  (ng)  and  the 

total  number  of  rays  -Ur,  and  pe  (size:  x  1)  is  the  sinogram  data  in  the  angle  0. 

Let  the  corrected  data  of  p  be  p'  which  consists  of 

=  Ag,r{Ng,r{Wr  •  Pe,r  -  R0,r)  -  Sg,r),  9  =  1,2,  ...  ,ng,  T  =  1 , 2,  .  .  .  ,  (B.l) 

The  terms  Ae,r  and  Ne,r  in  Eq.  (B.l)  are  multiplicative  correction  terms  for  atten¬ 
uation  and  detector  efficiency  (normalization).  Wr  corrects  for  the  fraction  of  time 
spent  in  each  wobble  position  and  is  generally  independent  of  angle.  Re,r  and  Sg,r 
are  the  estimates  of  random  and  scatter.  The  terms  pg,r,  p'g,r->  ^o,r,  and  Se,r  have 
units  of  counts,  while  all  other  terms  are  dimensionless. 

Therefore,  for  any  pixel  Qi  with  index  i  at  position  Xi,  yi,  the  pixel  value  calcu¬ 
lated  by  the  EBP  is  the  convolution  of  the  corrected  projection  data  Pe,r  with  the 
reconstruction  filter  /i,  that  is, 

ng  nr  . 

=  h{^i  COS  9  +  Hi  sin  0  -  (B.2) 


In  practice,  additional  linear  interpolation  steps  are  applied  and  the  reconstructed 
pixel  value  becomes 

ne  Ur 

*  =  (b-3) 

0=1  r=l 

where  is  unique  for  the  pixel  i  in  FBP  at  angle  6  and  ray  r. 

B.1.1  Variance  of  Pixel  in  an  ROI 


From  Eq.  (B.3),  the  variance  of  the  pixel  i  can  be  simply  written  as 


nr  .  .X  V  2 

Var{i?,}  =  £  £  •  Var{p'e,  J. 

0=1  r=l 


(B.4) 


In  a  PET  system,  it  is  usually  assumed  that  the  significant  source  of  noise  in  pg  ,.  is 
due  to  the  Poisson  statistics  in  the  raw  projection  data  pe,r,  that  is,  the  noise  in  the 
term  Ag^r,  Ne,r,  Wr,  Re,^  and  5e,r  is  small.  Also,  the  raw  count  pB,r  are  statistically 
independent  Poisson  distribution,  so  are  the  by  the  Eq.  (B.l).  Therefore,  the 
pixel  variance  in  Eq.  (B.4)  can  be  approximated  as 


TLq  Thr  r\ 

Varte}  ^  EEKO  •  •  Var(p,,) 

^=1  r=l 


(B.5) 


0=1  r=l 
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We  can  further  assume  that  for  each  angle  the  product  of  pg^r  ^.nd  \A0  ,.Nq,.W^\ 
is  relatively  uniform  in  the  ROI,  i.e.,  we  can  approximate  [A\^^Ng^r^r\P6,r  by  some 
representative  of  them,  namely  [A]N'^eW‘^]pe.  Therefore,  Eq.  (B.5)  becomes 


Varte)  = 

=  EKA'X^Jp.-EKT 

e=i  '•=1 

E  [aWw^]  p, 

0-1 

=  c,-'£[a\n1W‘]p, 


ne 

E 

0=1 


(B.6) 


where  ||h||2  is  defined  as  the  2-norm  of  the  filter  vector  h  and  he  is  the  filter  vector 
associated  with  the  angle  6.  We  used  the  observation  that  the  2-norm  of  the  filter 
vector  he  in  all  angle  are  approximately  equal  in  an  ROI,  that  is,  ||h|l2  «  llbelb?  for 
all  0.  Hence,  the  value  Cj  =  ||h||^  can  be  determined  a  priori  for  the  convolution 

filter  used  in  the  EBP. 
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B.1.2  Covariance  of  Pixels  in  an  ROI 

Following  the  similar  steps  in  deriving  the  variance  of  a  pixel,  the  covariance  of  any 
two  pixels  can  be  computed  as 


ne  Ur 


Cov{gj,gk)  =  E  E 

^=1  r=l 
ne  Ur 


(B.7) 


e=i  r=l 


where  gj  and  gk  are  two  pixel  values  of  interest  with  reconstruction  filters  and 
hfl,  respectively.  The  correlation  coefficient  of  pixels  gj  and  gk  is  given  by 


p{9ji9k^  ~ 


^ov{gj,gk) 


y'Var(firj)  •  Var(^fc) 


E  E 

e=l  r=l 


p,. 


J{Z  £  (41^)'  ■  { I,  I. 

(B.8) 


Assume  that  the  pixel  variances  of  and  gk  in  an  ROI  are  similar,  and  that  the  prod¬ 
uct  A|  r^erW^  •Pe,r  is  relatively  uniform  for  those  projection  lines  that  substantially 
contribute  to  the  sums,  Eq.  (B.8)  can  be  further  simplified  as 


pi9ji9k)  — 


2  £  4545 

e=l  r=l 


E  E  (45)1  [2  £  (45 

e=l  r=l  ^  '  J  L0=1  r=:l  ^ 


wy 
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(B.9) 


E  E  1^)1^ 

6=1  r=l _ 

TIQ  Tljn 

E  E/A.,r)^ 

9=1  r=l 


where  C2  =  l/(ci  •  ne).  Thus,  the  correlation  coefficient  of  any  two  pixels  in  an  ROI 
is  independent  of  the  raw  data  and  Eq.  (B.9)  can  be  used  for  a  set  of  pixel  pairs  at 
various  positions  to  produce  a  p{d)  table. 


B.2  Covariance  Computation  in  MAP  images 


The  method  used  in  this  dissertation  to  compute  the  covariance  of  the  MAP  re¬ 
construction  image  was  originally  developed  by  Fessler  [19,  20].  Recently,  Jinyi  and 
Leahy  [31]  proposed  a  fast  algorithm  to  reduce  the  computational  complexity  in 
Fessler’s  method.  The  following  is  a  brief  review  of  the  fast  covariance  computation 
of  MAP  image  described  in  [31]. 

B.2.1  MAP  Reconstruction 

A  PET  system  can  be  modeled  as  independent  Poisson  random  variables  y  with 
mean  y  related  to  the  emission  intensity  image  x  through  an  affine  transform 


y  =  Px  +  n 


(B.IO) 
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where  P  is  the  detection  probability  matrix  and  n  accounts  for  the  presence  of  noise 
in  the  data.  The  detection  probability  matrix  is  modeled  as  P  =  D[ni]G,  where 
D[ni]  is  a  diagonal  matrix  containing  the  correction  factors,  n,,  i.e.  the  product  of 
the  detector  normalization,  dead-time  and  attenuation  correction  factors^.  G  is  the 
geometric  projection  matrix  representing  the  probability  that  an  emission  from  each 
voxel  in  the  image  produces,  in  the  absence  of  attenuation  effects,  a  photon  pair  at 
each  of  the  detector  pairs  in  the  system.  The  log-likelihood  function  for  the  Poisson 
data  model  is 

L{y\^)  =  -  yi  +  iogyi!  (B.ii) 

t 

The  MAP  reconstruction  is  usually  described  as  the  maximizer  of  the  log  posterior 
probability: 

x{y)  =  B,TgmaxL{y\x)  -  (3U{x)  (B.12) 

where  /?  is  the  hyper-parameter  that  determines  the  relative  influence  of  the  prior 
and  likelihood  terms,  L{y\x)  is  the  log-likelihood  function  for  the  Poisson  data  model, 
and  U{x)  is  the  Gibbs  energy  function  depending  on  a  sparse  neighborhood  matrix 

C. 

^D{ni]  denotes  a  diagonal  matrix  with  the  (z,i)th  diagonal  element  equal  to  n,-. 
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B.2.2  Approximation  of  Covariance 

Based  on  a  fixed  point  of  the  MAP  objective  function  and  using  a  truncated  Taylor 
series  expansion  of  the  implicit  estimator,  Fessler  [19,  20]  developed  a  closed  form 
for  the  computation  of  covariance  in  the  MAP  reconstruction  image 

Cov(x)  ^  [F  +  /?C'C]-'F[F  +  ^C'C]-'  (B.13) 

where  F  =  P'I>[l/y,]P  is  the  Fisher  information  matrix.  This  closed  form  involves 
a  computation  of  the  inverse  of  a  Hessian  matrix  or  solving  a  related  set  of  linear 
equations.  Hence,  to  reduce  the  computation  complexity,  Jinyi  and  Leahy  [31]  pro¬ 
posed  a  fast  computation  method  which  uses  the  following  approximation  for  the 
Fisher  information  matrix 

F  ^  D«G'GD«  (B.14) 


where  =  D[Kj]  with 

Then, 

Cov(i)  “  d;‘b(i)-‘g'gb(x)-‘d;' 


(B.15) 


(B.16) 
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where 


B{x)  =  G'G  +  /3D:^aCB:\  (B.17) 

The  covariance  can  be  further  approximated  using  the  following  observations:  (i) 
the  matrix  B(a:)  represents  a  local  image  blurring  operator  since  the  prior  energy 
function  is  defined  on  a  local  neighborhood  and  the  blurring  function  G'G  is  also 
local;  (ii)  the  correlation  between  voxels  drops  off  rapidly  as  a  function  of  the  distance 
between  them.  The  covariance  with  respect  to  voxel  j  is  therefore  dominated  by  the 
contribution  of  kj  and  the  following  approximation  can  be  made: 

Cov,(x)  ^  D;'B(x)-'G'GB(x)-'D;'e,- 

^  KfKi^Kf)-^G'GK{/3Kf)-^ej  (B.18) 

where 

K{f3Kf)  =  G'G  +  pKfC'C.  (B.19) 

Since  G'G  and  C'C  correspond  to  the  shift  invariant  blurring  operators,  they  have  a 
block  Toeplitz  structure  and  can  be  approximately  diagonalized  using  the  2D  discrete 
Fourier  transform,  i.e., 

G'G  ^  Q'D[A,]Q 
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C'c  =  Q'D[fjii]Q 


(B.20) 


where  Q  and  Q',  respectively,  represent  the  Kronecker  form  of  the  2D  DFT  matrix 
and  its  inverse,  A^’s  are  the  2D  Fourier  transform  of  the  system  geometric  response 
at  the  central  voxel,  and  /li’s  are  the  2D  Fourier  transform  of  the  central  column  of 
C'C.  By  Eq.  (B.20),  Eq.  (B.19)  can  be  written  as 


Then,  the  final  covariance  approximation  formula  becomes 


Cov,(i)  ^ 


.-2 


Q'D[ 


A. 


(Ai  +  ^/it) 


rlQ^r 


(B.22) 


In  cases  where  only  a  small  ROI  is  considered,  the  correlation  structure  can 
be  assumed  to  be  invariant  within  the  region.  Hence  the  covariance  matrix  can  be 
computed  by  one  modified  backprojection  to  compute  the  /Cj’s  and  one  2D  FFT.  One 
advantage  of  the  above  derivation  is  that  A,-  and  fXi  are  independent  of  the  data  and 
can  be  pre-computed.  Another  advantage  of  Eq.  (B.22)  is  that  it  is  readily  inverted 
which  is  an  attractive  property  for  the  pre-whitening  process  of  the  non-white  PET 
data  for  computer  observer  lesion  detection. 
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Appendix  C 

GLRT  Can  Not  Be  Improved  by 
Orthogonal  Subspace  Refining 


Given  two  identified  lesion  subspace  (H)  and  normal  tissue  subspace  (S),  the  derived 


GLRT  in  chapter  5  can  be  written  as 


Ls{y) 

Lr{y) 


y'^^Hsy 

y^Pgy 

y^P^y 


(superimposed  model) 
(replacement  model). 


(C.l) 


The  refined  subspace  matrices  by  the  orthogonal  projection  of  (H)  onto  the  comple¬ 


ment  subspace  of  (S)  can  be  written  as 


S'  =  S,  H'  =  P^H 


(C.2) 
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where  represents  the  orthogonal  projection  onto  the  complement  subspace  of 
(S).  Note  that  S'  is  orthogonal  to  H'.  Then  we  have  the  following  properties: 


H'  =  P^H  =  H-PsH  Ph'  =  Ph  -Ppsn 

Ph>S'  =  Ph'  +  Ps'  =  Pp^,H  +  Ps'  =  Pp|H  +  Ps  =  PhS-  (C.4) 


Superimposed  model 

From  the  property  in  Eq.  (C.4),  the  GLRT  for  superimposed  model  can  be  written 
as 


^'(y) 


y^P^,y  _  y^PsY 
y^Ph's'Y  y^Pnsy 


^.(y) 


(C.5) 


wherethe  GLRTs  before  and  after  the  orthogonal  projection  of  the  lesion  subspace 
onto  the  orthogonal  complement  of  the  normal  tissue  subspace  are  shown  to  be 

equal. 

Replacement  model 

The  GLRT  for  the  replacement  model  using  H'  and  S'  is 


L'iy)  = 


y'^P^y 

y^Pn-y’ 


(C.6) 
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Then,  by  using  the  property  in  Eq.  (C.3),  we  have 


1  ^  y’'PH'y 

i;(y)  y^P^.y 

[Pff  +  Ppsh]  y 

"  y^'P^.y 

1  y^PpsHy 

“  i.(y)  y^'P^T  ■ 


(C.7) 


Hence, 


K(y)  = 


1  I  y^rp5H> 

Lr(y)  y'^P^.y 


<  Lr{y)- 


(C.8) 


Eq.  (C.8)  shows  that,  for  the  replacement  model,  the  performance  of  the  GLRT 


after  applying  the  orthogonal  projection  refined  subspaces  is  inferior  to  that  using 
the  original  subspaces,  because  the  scaler,  y^pf  y  ,  is  positive. 
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Appendix  D 

Multi-Pixel  GLRT  Derivation 


Given  an  ROI  with  P  number  of  pixels,  let  the  P-dimensional  column  vector  f,  = 
[fi{l),fi{2),  ■  ■  ■,fi{P)V  denote  the  spatial  pixel  data  for  the  i-th  frame  and  let  the 
ROI  data  Fpxiv  =  We  assume  that  all  the  columns  of  the  matrix  F 

are  independent.  Let  YnxP  =  [yi5y2)  •  •  •  ^yp]  =  ■ 

The  hypothesis  data  model  for  the  ROI  under  the  replacement  model  assumption 

is 

f  Hq  :  Y  =  Xo  +  No 

\^i:Y  =  Xi+Ni 

where  the  noise  N,-  is  assumed  to  be  a  multivariate  normal  (MVN)  with  zero  mean. 
It  is  also  assumed  that  the  signal  Xi  obeys  the  linear  subspace  model,  i.e.. 


Xo  =  S$o 

Xi  =  H$i  (D.2) 


no 


where  H  and  S  are  the  lesion  and  normal  tissue  subspaces,  respectively,  and  is 
the  corresponding  coefficient  matrix.  Hence,  the  mean  of  temporal-spatial  ROI  data 
becomes 


E{Y}  = 


S$o,  under  hypothesis  Hq 
under  hypothesis  Hi. 


(D.3) 


We  assume  that  the  spatial  covariance  matrix  for  each  frame  has  the  same  structure 


but  with  different  energy  level  [7].  That  means  the  spatial  covariance  matrix  for  the 


i-th  frame  can  be  written  as  ctjM,  i.e., 


cov  I,- 


(D.4) 


where  M  is  the  spatial  inter-pixel  correlation  matrix  for  all  frames.  Let  Z  =  M  , 


then  the  variance  matrix  for  the  whitened  data  Z  is 


var(Z)  = 


ai  <J2  ctn 

ai  (72  •  •  •  (7JV 

(Ti  (T2  •••  (TiV 


(D.5) 


J  PxN. 


Let  R  =  [ri,  r2.  •  •  • ,  rp]  =  Z^.  Then  the  covariance  matrix  for  each  column  of  R  is 


cov(r,)  = 


(7i  0 

0  <72 

0  0 


0 

0 

CTN 


^A,  ,  =  1,2, ■■■,/>.  (D.6) 


NxN 
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The  ROI  data  will  be  statistically  described  as  a  Gaussian  probability  density  funo 


tion  (PDF)  based  on  the  temporal-spatial  data  matrix  R.  After  the  whitening 
process,  by  the  independence  of  the  columns  r,  in  R,  the  joint  Gaussian  PDF  of  R 
under  hypothesis,  say  Hi,  becomes 


/i(R) 

=  /(ri,r2,---,rp|ifi) 

=  n^ifivilHi) 

=  nf,,(2^)-'"^|A.|-‘/^exp  [r,  -  E(r,)|^  (r.  -  E(ri)l} 

=  (2x)-'"’/^|Ax|-'’'^exp|^f:iri  -  E(rOrAr‘  [r,  -  E(r.)l| 

=  (2x)-'^'’/"|A.|-'’/"etr  [R-  E(R)f  Af*  |R  -  E(R)]} 

=  (2ir)-"'’'"|Air'’'"etr{^Ar‘[R-E(R)l[R-E(R)f}  (D.7) 

where  “etr”  denotes  “exp  trace”,  and 

E(R)  =  E  (Z^)  =  E  =  E 

=  E(Y)M-'/"  =  H^iM-^/1  (D-8) 
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If  we  further  assume  that  the  mean  of  image  is  the  same  for  all  frame,  i.e., 

(D-9) 

where  </)i  is  the  coefficient  vector  for  all  frames  and  1^  =  [1, 1,  •  •  • ,  1].  Then, 

E(R)  = 

^  bxs^  (D.IO) 

where  bi  =  and  Hence,  the  PDF  for  R  becomes 

/(R)  =  (27r)-'"’'*|A.r''/^etr|:^Ar'(R-b,sJ')(R-b,s’’)’'}.(D.ll) 

In  order  to  apply  the  generalized  likelihood  ratio  test  (GLRT)  for  the  matched 
subspace  detector,  the  next  step  is  to  find  the  maximum  likelihood  estimate  (MLE) 
of  Ai  and  cf)^.  Based  on  (D.ll),  the  MLE  of  Ai  can  be  simply  written  as 

A.  =  i(R-b.s’-)(R-b.sT-  (D12) 

The  MLE  of  (f>i  can  be  derived  from  the  minimization  of  |Fbi  |,  i-e., 

jL(R)  =  min|Fbi|  (D.13) 

bi 
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where 


Fbi  =  (R-biS^)(R-biS^)^.  (D.14) 

Now  expand  Fb,  as  follows: 

Fb,  =  RR^  -  Rsbf  -  bjs^R^  +  bis^sbf .  (D.15) 

Let  s  =  (s^s)i/'s  and  bi  =  (s^s)i/2bi,  so  biS^  =  b^s^.  After  the  change  of 
variables, 

Ff,  =  RR^  -  Rsbf  -  bis^R^  +  bibf 

bi  ^ 

=  (bi  -  Rs)(bi  -  Rs)^  +  RR^  -  (Rs)(Rs)^ 

=  (bi  -  Rs)(bi  -  Rsf  +  G  (D-16) 

where  G  =  RR^  -  (Rs)(Rs)^  and  G  can  be  shown  to  be  invertible  [52].  Hence, 

|F£,J  =  |G|-|I+(bi-Rs)(bi-RsfG-'| 

=  |Gl-|I  +  aa^G-'|  (0-17) 

where  a  =  bi  -  Rs.  Now  multiply  the  matrix  D  =  1  +  aa^G"^  by  a  as  follows: 

Da  =  a  +  a(a^G'^a) 
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=  (1  +  a^G  ^a)a 


(D.18) 


which  shows  that  a  is  an  eigenvector  of  D  with  the  associated  eigenvalue  1  +  a^G  ^a. 
Lgl;  yi,  for  i  —  1,2)  ’ ' '  f  JV  —  1,  be  some  orthogonal  set  of  vectors,  all  perpendicular 
to  the  vector  G"^a.  Then, 


Dvi  =  (I  +  aa^G"^)vi 

=  V, •  + a  [(G"^a)^v,] 

=  Vi 


(D.19) 


which  means  Vj  is  an  eigenvector  of  D  with  eigenvalue  equal  to  one,  for  i  —  1,2,  •  •  • ,  N 
1.  So  Eq.  (D.16)  becomes 

IFbJ  =  |Gl.(l  +  a^G-^a) 

=  |G|- [l  +  (bi-RsfG-'(bi-Rs)] 

=  |G|  •  [l  +  (H<^i  -  -  Rs)]  (D.20) 

where  H  =  (s^s)^/^H.  The  minimum  of  I  can  be  derived  by  using 

^m^Qm  =  2  ( Qm  (D.21) 
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provided  Q  is  independent  of  (py.  Then  we  have 


=  (H^G-'H)-‘H’'G-‘RI.  (D.22) 

Hence, 

|min 

=  |G|  •  |l  +  [H(H^G-^H)-'H^G-^Rs  -  Rs]^  C 

•  [h(H^G-'H)-'H^G-'Rs-Rs]} 

=  |G|  •  |l  +  [G-^/2H(HTG-iH)-iH^G-^/2G-i/2i^g  _  q-i/^Rs]"’ 

•  ■g-^/^H(H^G-^H)~^H^G“^/^G-^/^Rs  -  G-^/^Rs]} 

=  |G1 .  {l  +  [H(H^H)-^H^G-i/2Rs  -  G-^/2rs]^ 

•  [h(H^H)"^H^G"^/^Rs-G-^/^Rs]} 

where  H  = 

=  |G|  ■  {l  +  [(Pb  -  I)G-‘/^Ri]^  [(Pb  -  I)G-''^Rl]} 

=  |G|  ■  ll  +  [G"'/^Rs]’^Pfl  [G"‘^^Rs]J  (D.23) 

where  Pfl  is  an  orthogonal  projection  matrix  onto  the  subspace  of  H,  and  Pg  is  an 
orthogonal  projection  matrix  onto  the  complement  subspace  of  H. 


Now,  the  GLRT  for  the  ROI  data  matrix  can  be  derived  easily.  Based  on  the 
assumptions  made  for  the  ROI  data  model  in  this  Appendix,  the  detection  problem 
becomes  a  test  of  the  distributions; 


f  Ho 

:  Y:  N 

S<f)ol'^,Ao  ®  M 

:  Y:  N 

Ai  0  M 

(D.24) 


which  is  equivalent  to 


i?o  :  R  =  YM-^/2  .  jv 
Hi  :  R  =  YM-^/2  :  N 


(D.25) 


where  ®  denotes  the  Kronecker  product.  The  likelihood  ratio  test  can  be  written  as 


=  L{huAi-,R) 
L(<^o>Ao;Y)  I(bo,Ao;R)‘ 


(D.26) 


Hence,  the  GLRT  can  be  derived  by  substituting  the  MLE  of  A,  and 

=  .  (D.27) 

^  ’  L(bo,Ao;R)  VlAii; 

Then  it  is  more  convenient  to  replace  the  GLRT  by  the  (P/2)-root  GLRT 

|G|  •  {l  +  (G-i/2Rg)Tp±(G-i/2Rg)J. 

|G|  •  {l  +  (G-i/2Rs)^P^(G-i/2Rs)} 

l  +  (G-^/WPi(G-^/^Rs)  281 

-  1  +  (G-i/2Rs)^P^(G-i/2Rs)‘  ^  ^ 
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Appendix  B:  Presentation  at  SNM  Annual  Conference, 
Torrance,  Ca.  June  2001 
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Assessment  of  ROI-based  Time 
Activity  Anaiyses  in  Dynamic 
PET  for  Oncoiogy 

X.  Yu,  C.C.  Huang  and  P.  S.  Conti 


OBJECTiVE 

•  To  develop  a  method  for  detecting  tumors 
in  PET  images  based  on 

—  time  activity  curve  (TAC)  information 

•  To  perform  an  assessment  on  three  TAC 
evaluation  methods  for  tumor  detection 
with 

—  receiver  operating  characteristic 
(ROC)  study 


APPROACH 


1)  Implement  three  TAC  evaluation 
methods  in  computer  software 

2)  Build  a  machine  observer  with  three 
different  TAC  evaluations 

3)  Generate  phantom  dynamic  raw  data 
with  or  without  artificial  lesions  and 
reconstruct  FBP  images  from  the  raw 
data 

4)  Use  machine  observer  to  evaluate  the 
phantom  images  and  perform  ROC  study 


TAC  EVALUATION  METHOD  1 

•  Assuming  no  spatial  correlation 
between  pixels 

—  Averaging  TACs  over  ROI 

—  Averaging  operation  will  lower  the 
noise  to  certain  extent  if  no  inter- 
pixel  correlation 


TAC  EVALUATION  METHOD  2 


•  Assuming  inter-pixel  correlation 

—  Compute  ROI  covariance  matrix 
frame  by  frame 

—  Averaging  TACs  over  ROI  after 
inter-pixel  decorrelation 


TAC  EVALUATION  METHOD  3 


•  Assuming  that  each  frame  has  the  same 
inter-pixel  correlation  structure  R,  but  with 
different  energy  levels:  Rj  =  Oj  R 

—  Compute  ROI  covariance  matrix  from 
the  least  noisy  frame 

—  Decorrelate  ROI  with  the  covariance 

—  Averaging  TACs  over  ROI  after 
inter-pixel  decorrelation 


COVARIANCE  MATRIX 
COMPUTATION 

Huesman’s  analytical  computation: 

—  based  on  uncorrelation  of  PET 
projection  data  and  linearity  of  FBP 

Carson’s  simplification  formulae: 

—  assume  that  the  statistics  of 
neighboring  projections  are  similar 

—  covariance  is  independent  of  raw 
data 


MACHINE  OBSERVER 

Design  principle:  the  maximum 
likelihood  ratio  test 

Function:  For  a  given  ROI,  make  a 
decision  based  on  the  likelihood  that 
the  observed  TACs  are  in 

-  a  lesion  or 

-  normal  tissues 


DATA  FOR  STUDAY 


Computer  generated  phantom 
dynamic  data 

-  Each  data  sequence  was  set  to  contain 
5  artificial  lesions  or  no  lesions 

-  50  sets  in  total  are  generated  for  each 
setting 


PHANTOM 


PHANTOM:  FBP 


IMAGE  GENERATION  AND 
RECONSTRUCTION 

Forward  projection  with 
“  system  blurring 
“  Poison  noise  generation 
Filtered  backprojection  (FBP)  with 

-  Han  filter  and  0.5  cutting 
frequency 


PHANTOM:  FBP 


NOTATION 


TP: 

True  Positive 

TN: 

True  Negative 

FP: 

False  Positive 

FN: 

False  Negative 

PPV: 

Positive  Predication  Value 

NPV: 

Negative  Predication  Value 

DEFINITIONS 


Accuracy: 

Sensitivity: 

Specificity: 

PPV: 

NPV: 

Az: 


(TP+TN)/(TP+TN+FP+FN) 

TP/(TP+FN) 

TN/(TN+FP) 

TP/(TP+FP) 

TN/(TN+FN) 

Area  under  curve 


ROC  PERFORMANCE 


Method  1 

Method  2 

Method  3 

Accuracy 

0.79 

0.77 

0.88 

Sensitivity 

0.77 

0.77 

0.79 

Specificity 

0.81 

0.77 

0.97 

PPV 

0.80 

0.77 

0.96 

NPV 

0.81 

0.77 

0.97 

Az 

0.93 

0.92 

0.98 

RESULTS  AND  COMPARISON 


•  Compared  to  Methods  1  and  2,  Method 
3  increased: 

—  Specificity  by  20% 

—  Accuracy  by  1 1  % 

—  Area  under  curve  by  6.5% 
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CONCLUSION 

•  No  significant  difference  between  Methods  1 
and  2,  possibly  due  to 

—  inaccurate  covariance  matrix 
estimation  in  the  early  noisy  frames 

—  simplified  covariance  matrix  computation 

•  Constrained  inter-pixel  decorrelation 
improves  TAC-based  lesion  detection 


